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We study the thermodynamic properties of the three-dimensional Blume-Capel model on the 
simple cubic lattice by means of computer simulations. In particular, we implement a parallelized 
variant of the multicanonical approach and perform simulations by keeping a constant temperature 
and crossing the phase boundary along the crystal-field axis. We obtain numerical data for several 
temperatures in both the first- and second-order regime of the model. Finite-size scaling analyses 
provide us with transition points and the dimensional scaling behavior in the numerically demanding 
first-order regime, as well as a clear verification of the expected Ising universality in the respective 
second-order regime. Finally, we discuss the scaling behavior in the vicinity of the tricritical point. 

PACS numbers: 75.10.Nr, 05.50.-|-q, 64.60.Cn, 75.10.Hk 


I. INTRODUCTION 

The Blume-Capel (BC) model consisting of a spin-one 
Ising Hamiltonian with a single-ion uniaxial crystal field 
anisotropy Hi is one of the most studied models in 
the communities of Statistical Mechanics and Condensed 
Matter Physics. This is not only because of the rela¬ 
tive simplicity with which approximate calculations for 
this model can be carried out and tested, as well as the 
fundamental theoretical interest arising from the richness 
of its phase diagram, but also because versions and ex¬ 
tensions of the model can be applied for the description 
of many different physical systems, some of them being 
multi-component fluids, ternary alloys, and ^He - ^He 
mixtures (3|. Recent applications of the BC model in¬ 
clude analyses of ferrimagnets, as discussed in a thorough 
contribution by Selke and Oitmaa Q. 

The BC model is described by the Hamiltonian 

H = - J ^ = Ej + AEa, (1) 

(b) » 

where the spin variables cri take on the values { — 1, 0, -1-1}, 
{ij) indicates summation over all nearest-neighbor pairs 
of sites, and J > 0 is the ferromagnetic exchange interac¬ 
tion (here we set J = 1 and fee = 1 to fix the temperature 
scale). The parameter A is known as the crystal-field 
coupling that controls the density of vacancies (ui = 0). 
For A —)■ —oo vacancies are suppressed and the model 
maps onto the Ising model. We always employ periodic 
boundary conditions. Note here that the second formu¬ 
lation of the Hamiltonian via the definitions of Ej 
and Ea , will allow us to define the necessary observables 
for the application of our finite-size scaling (FSS) scheme 
that will be discussed in detail below. 

As is well known, the pure and disordered ver¬ 
sions of the model of Eq. © have been analyzed, be¬ 
sides the original mean-field theory 0, i, by a va¬ 
riety of approximations and numerical approaches, in 
both two (2D) and three dimensions (3D). These in¬ 
clude the real-space renormalization group Q, Monte 



FIG. 1: General sketch of the phase diagram of the 3D Blume- 
Gapel model in the temperature - crystal field plane. The 
phase boundary separates the ferromagnetic (F) from the 
paramagnetic (P) phase, in which the solid line indicates con¬ 
tinuous and the dotted line first-order phase transitions, re¬ 
spectively. The two lines merge at the tricritical point (TP), 
as highlighted by the black rhombus. The limiting cases of 
T = 0 and A = 0 are marked on the relevant axis with Ao = 3 
and To = 3.195(1) [l^l, respectively. The horizontal arrows il¬ 
lustrate the direction of crossing the phase boundary at fixed 
temperatures, studied in this work, in both the second-order 
(Ti = 2.0) and first-order (T 2 = 1.0 and T 3 = 0.9) regimes, as 
well as in the vicinity of the tricritical point (Tt = 1.4182) [l4| . 


Carlo renormalization-group calculations Q , e-expansion 
renormalization groups Q, high- and low-temperature 
series calculations a phe nomenological FSS analysis 
using a strip geometry [^, [l^ , and of course Monte Carlo 
simulations [lll - l^ . 

The phase diagram of the pure model consists of a 
segment of continuous Ising-like transitions at high tem¬ 
peratures and low values of the crystal field. This ends 
at a tricritical point (At,Tt) where it is joined with 
a second segment of first-order transitions ending for 
T = 0 at Aq = zJ 14, with z denoting the coordina¬ 
tion number of the considered lattice. In the present 
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case of a simple cubic lattice, where z = 6 , it follows 
that Aq = 3. A general sketch of the phase diagram 
is given in Fig. [T] and is outlined in the relevant cap¬ 
tion. The location of the tricritical point, marked by 
the black rhombus in Fig. [U has been estimated by De- 
serno 0 using a microcanonical Monte Carlo approach 
to be (At,Tt) = (2.84479(30), 1.4182(55)). 

The scope of the present paper is to present a com¬ 
plementary study of the 3D BC model embedded in 
the simple cubic lattice. Our simulations follow a so¬ 
phisticated numerical scheme, outlined in the following 
Sec.in using as a platform the multicanonical approach. 
This is especially suitable for the study of systems that 
undergo a first-order phase transition, where it is well- 
known that numerical simulation is a hard task to per¬ 
form. One interesting aspect of the present work is that 
we cross, in our simulations, the phase boundary of the 
system along the crystal-field axis, keeping the tempera¬ 
ture fixed. Thus we obtain relevant thermodynamic ob¬ 
servables as a function of the crystal field A and we per¬ 
form a FSS analysis on a different basis. This analysis is 
presented in Sec. imi in both the first- and second-order 
regimes of the model, where the second-order regime is 
used as a test case of our scheme and more attention 
is paid to the first-order regime of the phase diagram, 
which is computationally much more challenging. In par¬ 
ticular, we consider three different temperatures, one in 
the second-order regime, Ti = 2 . 0 , and two of them in 
the first-order regime, T 2 = 1.0 and T 3 = 0.9. More¬ 
over, we discuss the scaling properties in the vicinity 
of the proposed tricritical point, by performing addi¬ 
tional simulations and relevant analysis at the temper¬ 
ature Tt = 1.4182 suggested in Ref. [3. This contribu¬ 
tion is ended in Sec. IIVI where a brief summary of our 
conclusions is given together with an outlook for future 
work. 


II. NUMERICAL METHOD AND SCALING 
OBSERVABLES 

We apply a multicanonical method [ 2 ^, with the 
slight modification to yield a flat histogram not in the 
total energy A, but rather in E^. The multicanonical 
method allows one to increase the probability to sam¬ 
ple otherwise suppressed states and, with it, overcome 
emerging barriers. Hence, it is an optimal tool to study 
first-order phase transitions. The canonical expectation 
value weights all observables of the phase space with the 
Boltzmann weight, which leads to the general form 

= ( 2 ) 

X 

where Z = jg partition sum, /3 = I/T is 

the inverse temperature, and x stands short for the spin 
configurations. For the usual multicanonical method, 
the Boltzmann weight in the canonical probability dis¬ 


tribution exp{—f3E(x)} is replaced by a weight func¬ 
tion W{E{x)), that is iteratively modified to yield a 
flat energy histogram. At this point, we can rewrite 
E = Ej + AAa, separate the probability distribution, 
and replace the Boltzmann weight depending on Ea- 

^-PEj g-/3AE^ ^ ^-pEj ^ ^ (3^ 

Considering a fixed inverse temperature /3, one is then 
able to iteratively adapt W (E/^) in order to yield a flat 
histogram in E/\. This is in fact quite similar to mul- 
timagnetic simulations and also suited for the applica¬ 
tion of a parallel implementation of the multicanonical 
method . We made use of this parallelization with up 
to 36 cores, which speeds up the iteration process and 
provides 36 independent production runs. The canonical 
expectation values at a certain point (/3, A) may then 
be estimated with standard histogram and time-series 
reweighting techniques [ 2 ^. Since the multicanonical 
simulation is still an importance sampling Markov chain, 
one only needs to consider the multicanonical variable 
Ea, illustrated for the case of time-series reweighting of 
a given observable O: 

' E.e-/5^^'^(DlV-i(^^A(x)) ■ 

Here, (• ■ •) clearly refers to the estimator of the expec¬ 
tation value and we will drop the subscripts in the fol¬ 
lowing. Figure [2] shows an example of the reweighted 
probability distributions of ca = £^a/H, where V = 
and L denotes the linear system size, at the transition 
field Aeqh, i-C. where the distribution shows two peaks of 
equal height. Well inside the first-order regime the sys¬ 
tem shows a barrier increasing with the system size, char¬ 
acteristic of the nature of the transition, which reaches 
at T 3 = 0.9 already for a system size oiV = 24^ spins the 
order of 10“^^° , see Fig. [2] (a). On the other hand, at the 
proposed tricritical point (see panel (b) with numerical 
data obtained at T = Ft = 1.4182 [T^) no definite judge¬ 
ment can be made. We observe that the distributions 
still show a double-peaked structure, yet with a much 
smaller barrier, which does not diverge with increasing 
system size. A more illuminating discussion of this spe¬ 
cial temperature is given below at the end of Sec. IHII 

It is worth noting here that the E^ are always in¬ 
teger in the range [0,U], which is a major advantage 
for the application of the present method. This would 
not be the case for the usual multicanonical [H, [2^ and 
Wang-Landau methods at fixed A based on the to¬ 
tal energy E, because E G M. for non-integer values of 
A. Additionally, our previous experience with the appli¬ 
cation of the Wang-Landau method to the BC model in 
the low-temperature regime for a fixed value of A sug¬ 
gested that one needs to be extremely careful in the de¬ 
tails of implementation. For instance, the need to simu¬ 
late large enough system sizes leads to the inevitable ap¬ 
plication of a multi-range approach of the Wang-Landau 
method, that is splitting the energy range in many subin¬ 
tervals [13 ■ This approach, although appearing to be 
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FIG. 2: (color online) Illustration of the reweighted probability distribntion at the transition field Aeqh with respect to 
BA = EaIV at (a) T = 0.9, that is well inside the first-order regime of the system (note the strongly increasing barrier 
with the system size), exhibiting a major suppression of intermediate states which is common for first-order transitions, and 
(b) at T = Tt = 1.4182 which is in the vicinity of the tricritical point. 


much faster than the straightforward one-range imple¬ 
mentation, may give rise to several problems with respect 
to the breaking of ergodicity of the process [28l - [30l| and 
possible distortions (systematic errors) induced on the 
density of states 0 . On the other hand, the numerical 
framework developed and applied in the current paper 
does not suffer from this type of inherent problems. On 
the contrary the parallelized variant of the multicanon- 
ical approach used, combined with the orthogonal scal¬ 
ing of the phase boundary has proven to be a promising 
scheme for the first-order transition regime. This will be 
clearly shown in the following section with the accurate 
estimation of transition points A* even for temperatures 
T < 1, which is a, commonly accepted, harsh numeri¬ 
cal task. Still, for the second-order regime both the cur¬ 
rent approach and any other type of generalized ensemble 
sampling method would give comparable results within 
the statistical errors as we have already verified by our 
preliminary numerical tests. In fact, the modification to 
a flat-histogram method in a sub-energy is not restricted 
to the present model or a spin system in general and has 
been applied in a similar way also to a polymer system in 
disorder [s^. Moreover, the formulation for other gener¬ 
alized ensemble methods is straight forward. 

In order to obtain transition points for the FSS anal¬ 
ysis, one usually considers the peak of the specific heat 
C, magnetic susceptibility X: or any other suitable tem¬ 
perature derivative of an order parameter . In princi¬ 
ple, the magnetic properties of the system show a more 
reliable behavior when one is interested in obtaining ac¬ 
curate estimates of critical points and it is a common 
practise along these lines to firstly estimate the magnetic 
exponents /3, 7 , and u, and then via the hyperscaling re¬ 
lation the exponent a of the specific heat. However, also 
other scaling approaches based on a different philosophy 
have been successfully used in the literature, depending 
always on the direction of intersecting the phase bound¬ 
ary of the model under study. For instance, it has been 


shown that for the 3D random-field Ising model at zero 
temperature, the field derivative of the bond energy Ej 
defines a specific-heat-like quantity from which one may 
produce accurate estimates of the critical exponent ra¬ 
tio ajv and whose shift behavior defines correct critical 
points in the field - temperature plane [s^. 

For the present study, where we are crossing the phase 
diagram at fixed temperature along the crystal-field axis, 
we may as well consider instead of the standard defini¬ 
tions, the field derivative of the form 9/3A. The deriva¬ 
tive of the expectation value ([ 2 ]) yields 

® = -/3[(Oi?A) - (0)(i5A)] + , (5) 

which is similar to any specific-heat-like quantity because 
in general the observable is independent on the variable 
and the last term drops out. This is true for either Ei, 
where i = J or A, however, the total energy E = Ej + 
AEa is no longer independent on the field which leads 
to 


d{E) 

dA 


-^[{EEa) - {E){Ea)] +{Ea) 


d{Ej) 

dA 


+ A 


d{EA) 

dA 


+ (Ea). 


( 6 ) 


However, the last line suggests - expecting a critical or 
diverging behavior - that we may limit our consideration 
to either of the energy contributions. In fact, we consider 
here only the field derivative of the spin-spin interaction 
term 


C(A) = = -P[{EjEa) - {Ej){Ea)] . (7) 


Similar considerations may apply also for other suitably 
defined thermodynamic functions that could provide us, 
for instance, with estimates of magnetic exponents men¬ 
tioned above. Yet, this task goes beyond the scope of the 
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FIG. 3: Specific-heat curves as a function of the crystal field 
A for T — 2.0 and several system sizes in the second-order 
regime. Smooth curves typical of continuous transitions with 
a clear shift behavior are observed. 

present work where we focus on the first-order transition 
regime of the BC modei and target oniy at a quaiita- 
tive comparison to the expected Ising criticaiity in the 
second-order transition regime. In fact, this comparison 
becomes even more ciear via the use of the straightfor- 
wardiy defined specific-heat-iike quantity ©• 

Reievant piots of C'(A) can be found in the next section 
in Figs. [3] and [5] and wiii be discussed there. However, it is 
obvious from these iiiustrations that there exists cieariy a 
maximum of C’(A) that moreover shows a shift behavior 
as weii. Let us define now AJ^ as the crystai fieid vaiue 
at which C'(A) attains its maximum, and as we shall see 
below this defines a suitable pseudocritical, or pseudo¬ 
transition, parameter that carries in itself the approach 
to the thermodynamic limit in the second- and first-order 
regime of the model, respectively. Similarly, we denote 
by (72 the value of the specific heat at this pseudocritical 
point (7(A2). The value of A^ is numerically determined 
by calculating the second derivative of (Ej) with respect 
to A, analogous to the above calculations, and finding the 
zero crossing. For this, we apply a bisection algorithm 
with time-series reweighting for the full data set as well 
as w subsets, each excluding 1/w measurements, for the 
jackknife error calculation |34l |. The jackknife method is 
similarly then used for (7^. 

To sum up, using the above scheme we have performed 
simulations for the three temperatures shown in Fig. [T] 
and outlined in the introductiom as well as for the tricrit- 
ical temperature Tt = 1.4182 [IJ. In all cases, we consid¬ 
ered various linear sizes within the range L = 8 — 28. In 
principle, our numerical approach could as well simulate 
even larger system sizes, especially in the second-order 
regime. Still we have found it useful to optimize our 
code following the needs of a careful inspection of the 
first-order regime, for which linear sizes of the order of 
L = 28 are already quite large, taking into account that 
we are well into the low-temperature part of the phase 
diagram of the model. 


FIG. 4: Shift behavior ([8| of the pseudocritical fields A2 
at T = 2.0 obtained from the peak location of the specific- 
heat curves shown in Fig. [S] The solid line is a fit for linear 
sizes L > Lmin = 18, giving a stable, under Lmin changes, 
critical field value Ac = 2.523(6). The inset illustrates the 
extrapolation of the effective exponent v, obtained by varying 
the cutoff Lmin of the fits, to the thermodynamic limit. 


III. FINITE-SIZE SCALING ANALYSIS 

In this section, we present the main results of our con¬ 
tribution based on a FSS analysis of the numerical data 
obtained with the method outlined above. As a first step 
we shall consider the scaling in the second-order regime 
of the model (T > Tt) and in particular at the tempera¬ 
ture Ti = 2.0. Let us point out here that when it comes 
to the second-order transition regime of the BC model, 
the modified multicanonical method as implemented here 
is by no means the method of choice if one wants to 
obtain high-accuracy estimates of critical exponents (or 
critical points) and cannot compete against other cluster- 
update methods or more involved generalized ensem¬ 
ble schemes [s^, IsJ especially tailored to this situation. 
However, a qualitative study at this regime allows a di¬ 
rect comparison to the extensive and precise literature 
of the simple 3D Ising model, thus serving as a clear 
cut test of the proposed scheme. Additionally it justihes 
the results and conclusions drawn from our study at the 
first-order regime, presented later in this section for two 
temperature values (T 2 = 1.0 and T 3 = 0.9 in Fig. [J), 
which is not so easy to control given the huge energy 
barriers illustrated in Fig. ISIa). Finally, in the last part 
of this section, we will discuss the scaling behavior of our 
observables at an estimate of the tricritical point (l^ . 

Figure [3] displays the specific-heat-like curves defined 
in Eq. © as a function of the crystal field A for T = 2.0. 
Several system sizes up to L = 28 are shown which ex¬ 
hibit a clear shift behavior. This is further quantified in 
Fig. 01 which presents the FSS behavior of the pseudo¬ 
critical fields A 2 estimated as the locations where the 
specific heat attains its maximum. As usual, for second- 
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FIG. 5: Specific-heat curves as a function of the crystal field 
A at r = 1.0 and two system sizes L — 8 (dashed line) and 
12 (solid line). The inset illustrates a typical first-order-like 
specific-heat curve for a system with linear size L = 28 on a 
double logarithmic scale. 

order phase transitions a scaling behavior of the form 

A2 = A, + (8) 

is used in order to describe the approach to the thermo¬ 
dynamic limit. It appears that this method of extracting 
pseudo-critical points from the maxima of some properly 
defined thermodynamic quantity is capable of produc¬ 
ing accurate estimates for both the critical crystal field 
Ac and also the correlation-length exponent v, assum¬ 
ing that its behavior follows the observed shift behavior 
of our pseudocritical fields A^. It is well known from 
the general scaling theory that, even for simple models, 
the equality between the correlation-length exponent and 
the shift exponent is not a necessary consequence of scal¬ 
ing [s^. Of course, it is a general practice to assume 
that the correlation-length behavior can be deduced by 
the shift of appropriate thermodynamic functions. 

In fact, the solid line in the main panel of Fig. |4] rep¬ 
resents a fitting of the form (O, using as a lower cutoff 
the linear size Fmin = 18. We have performed this type 
of analysis for several values of Fmin within the range 
8—18 and keeping of course the upper system size fixed 
at Fmax = 28. For each of these fits we have estimated an 
effective value of the correlation-length’s exponent which 
is plotted in the inset of Fig.|4]as a function of the inverse 
lower cutoff, i.e., the parameter I/Lmin- A linear extrap¬ 
olation to the infinite-limit size provides an estimate of 
iz = 0.62(4), which within error bars is compatible with 
the Ising universality exponent i/ = 0.6304(13) [i^, as 
expected. Regarding the value of the critical field we 
obtain the estimate Ac(ri = 2.0) = 2.523(6), that re¬ 
mained quite stable under the switching of the lower cut¬ 
off during the fitting procedure. Thus, up to this point 
we have verified through a rather different, “orthogonal” 
route the expected Ising universality in the second-order 
phase transition regime of the 3D BC model. 



I /V 


FIG. 6: Shift behavior ([9| of the pseudocritical fields AJ, 
obtained from the peak location of the curves (a sample of 
which is shown in Fig. El) for both temperatures T = 1.0 and 
T = 0.9 in the first-order regime of the phase diagram. The 
obtained transition fields A* are given in the main text. The 
inset is a mere enlargement of the ~ L~‘^ approach of the 
pseudotransition points A^ to L —>■ oo for the temperature 
T = 1.0. 

We now move on to the main objective of this work, the 
discussion of the characteristics of the transition in the 
first-order regime and its dimensional scaling behavior. 
Let us point out here before discussing our findings that 
crossing the boundary at the first-order transition regime 
at a fixed temperature is an orthogonal approach to the 
fixed-field ansatz. One advantage in the case of the BC 
model is a broad temperature range with a first-order 
transition in comparison to a small A-range. 

As we already discussed above, we have obtained nu¬ 
merical data at two temperatures in the first-order regime 
of the model, T 2 = 1.0 and T 3 = 0.9. A nice illustration 
of the first-order character of the transition at these tem¬ 
peratures is shown in Fig. [5l where we plot the specific- 
heat curves obtained from Eq. ([7]) for T = 1.0 and two 
system sizes, L = 8 and L = 12 (and L = 28 in the in¬ 
set at a double logarithmic scale). Clearly, a sharp peak 
is observed which becomes much more pronounced with 
increasing system size. 

Following a similar analysis as above, we study now the 
scaling of the pseudocritical fields A^ obtained from the 
sharp specific-heat peaks. In this case we would expect 
a scaling of the form 

Al= A* +bL-'^, (9) 

where d = 3 the dimensionality of the lattice and A* the 
transition field. The above shift behavior of the pseu¬ 
docritical fields A*j^ for both temperatures T = 1.0 and 
T = 0.9 in the first-order regime of the phase diagram is 
shown in Fig. [5] as a function of the inverse of the volume 
of the system. The solid lines are linear extrapolations to 
the infinite-limit size (for a clearer illustration of the lin¬ 
ear behavior see the corresponding inset). The obtained 
transition fields A* are 2.944(5) and 2.964(6) for T = 1.0 
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FIG. 7: Simultaneous fitting of the specific-heat maxima at 
both temperatures in the first-order regime. The expected 
~ L'^ scaling behavior is obtained as can be clearly seen. 


and 0.9, respectively. 

As another important aspect of the first-order regime 
in the phase diagram of the model, we study the scaling 
of the specific-heat maxima in Fig. [T] In particular, we 
plot the FSS behavior of the peaks for both temperatures 
considered in this regime, where the two solid lines show 
a simultaneous fitting attempt of the form 
simultaneous meaning that the two data sets share the 
same exponent during the fitting procedure. Of course, 
in a standard first-order phase transition, the exponent 
p is expected to be equal to the dimensionality d of the 
system, that is 3 in our case. The result for the exponent 
p of a simultaneous fit to the data for both temperatures 
with x^/dof w 0.8 is p = 3.00(2), which is in excellent 
agreement with the theoretical expectation p = d = 3. 

Further to the above successful study of criticality in 
the second-order regime of the model, it is now clear that 
the numerical method and scaling approach implemented 
in the present paper is able to capture as well the first- 
order characteristics of the transition within a good ac¬ 
curacy. This latter fact is of particular importance as we 
are dealing with the low-temperature first-order regime 
of the BC model, where it is common knowledge that 
most numerical methods fail to produce reliable estimates 
of transition points and criticality. Thus, the current 
method could be easily stretched to produce an accurate 
approximation of the phase boundary line for values of 
the crystal field A within the regime [At, 3]. 

The multicanonical method allows us to directly esti¬ 
mate the barrier associated with the suppression of states 
during the first-order phase transition, as shown in Fig. [21 
Considering distributions with two peaks of equal height, 
i.e., two equally probably states, leads to the formulation 
of a free-energy like barrier in the i?A-space, 


B = 


2/3A 


In 


eqh 


( 10 ) 


where Pmax and Pmin are the maximum and the local 




L-2 


(N 




L-2 


FIG. 8: Scaling of the barrier height B for all considered 
temperatures T <Tt. This barrier may be associated with the 
transition from a spin-0 dominated to a spin-±l dominated 
regime (see Fig. [2]). Then B/L?' plays the role of an interface 
tension for spin-0 strips. 


minimum of the distribution P{E/^), respectively. The 
resulting barrier connects a spin -0 dominated regime (Aa 
small) and a spin-±l dominated regime {Ea large). This 
shows large similarities to the Ising (lattice gas) model 
and the according droplet/strip transitions [H, l40l| . 
Thus, the association with condensation and strip for¬ 
mation of spin -0 clusters seems natural and we would 
expect a scaling behavior in three dimensions as B jLp' = 
a -h ciL~‘^ -|- O possibly with higher-order correc¬ 

tions [4l|. Figure IH] shows B/L"^ as a function of L~^ 
for T 3 = 0.9, T 2 = 1.0 and Tt = 1.4182 with fits of the 














7 



5 10 15 20 25 30 


L 


FIG. 9: Critical aspects of the 3D Blume-Capel model at 
the tricritical point proposed in Ref. [TJ: T = Tt = 1.4182. 

(a) Shift behavior of AJ, obtained from the location of the 
specific-heat peaks. The solid line is a power-law fit of the 
form ([8ll for L > Lmin = 20. The inset illustrates the infinite- 
volume extrapolation of the correlation-length effective ex¬ 
ponent by varying the lower cutoff Lmin during the fittings. 

(b) Scaling behavior of the specific-heat peaks again for the 
larger system sizes. The corresponding inset shows an infinite- 
volume extrapolation of the effective exponent ratio a/v using 
the same procedure as in (a). 


data including the first correction and up to the third 
corrections. While higher-order corrections describe the 
systematic dependence of the data better, the L —>■ oo ex¬ 
trapolations are consistent within error bars for both fits 
yielding the estimates = 0.0774(1) and <72 = 0.0540(2) 
for T 3 = 0.9 and T 2 = 1.0, respectively. In the vicinity 
of the tricritial point, at Tt = 1.4182, the extrapolation 
yields cr Ri 0 indicating as expected that the interface 
tension vanishes in the thermodynamic limit. 

In the last part of this section we discuss some scal¬ 
ing results in the vicinity of the tricritical point of the 
3D BC model 0 . We have performed additional simula¬ 
tions by hxing the temperature at the tricritical estimate 
T = Tt = 1.4182, as suggested by Deserno 0 , cross¬ 
ing again the phase diagram along the crystal-field axis. 
The results and relevant FSS analysis are given in Fig. [HI 
where one can clearly observe the departure from the 
Ising second-order universality class to the tricritical one, 
at least in terms of the estimated critical exponents. In 
particular, in panel (a) of this figure we present the shift 
behavior of obtained from the location of the specific- 
heat peaks at the above defined temperature. The solid 
line is a power-law htting of the form (|5]) for T > Lmin = 
20 and the estimate we obtain for the relevant (tricrit¬ 
ical) crystal-field value is At = 2.8446(3). This latter 


value compares very well to the value 2.84479(30) pro¬ 
posed by Deserno, using an empirical scaling of the coor¬ 
dinates of a latent-heat-like quantity of the model. The 
inset of panel (a) illustrates correspondingly the infinite- 
limit size extrapolation of the correlation-length’s effec¬ 
tive exponent by varying the cutoff Tmin during the fit¬ 
ting procedure, as also performed in the analysis within 
the second-order regime of the model (see Fig. S]). The 
obtained value of v, that isv = 0.517(37), is clearly differ¬ 
ent to that of the standard second-order Ising universality 
class, and within error bars, compatible to the theoreti¬ 
cal expectation of the Ising tricritical universality value 
of zz = 0.5 l3l.[ig . This result indicates that the estimate 
of Deserno [ij for the location of the tricritical point 
in the temperature - crystal-held plane is indeed quite 
accurate, and secondly it provides a strong test in favor 
of the implemented numerical and scaling scheme of the 
present paper. Further to these results, we present in 
Fig. O (b) the scaling behavior of the specihc-heat peaks 
C’l, following the scaling law C£ . The solid line 

is a power-law ht of this form, again for the larger system 
sizes, and the corresponding inset illustrates the inhnite- 
volume extrapolation of the effective exponent ratio ajv. 
This analysis leads to an estimate ajv = 0.98(4), again 
very close to the expected Ising tricritical universality 
value of ajv = 1 [1,0. 


IV. CONCLUSIONS AND OUTLOOK 

In this manuscript, we have presented a numerical 
study of the three-dimensional Blume-Capel model de- 
hned on a simple cubic lattice. By implementing a vari¬ 
ant of the multicanonical method, we have performed 
simulations of the model keeping a constant tempera¬ 
ture and crossing the phase boundary along the crystal- 
field axis. In this way we have obtained numerical data 
for several temperatures in both the first- and second- 
order regime of the model, as well as in the vicinity of 
the tricritical point. A standard hnite-size scaling analy¬ 
sis, mainly based on a properly defined specific-heat-like 
quantity, provided us with precise estimates for the tran¬ 
sition points in both regimes of the phase diagram and 
with a clear verification of the expected ~ scaling be¬ 
havior and the Ising universality class in the first- and 
second-order regimes of the model, respectively. 

An interesting feature of our study is related to the 
fact that we have been able to probe efficiently the low- 
temperature first-order regime of the phase diagram of 
the Blume-Capel model, a rather tricky numerical task, 
and obtain accurate estimates of transition points in the 
regime of strong crystal fields. Using the multicanonical 
method, and hence simulating otherwise strongly sup¬ 
pressed states, allowed us to measure the associated free- 
energy like barrier in the first-order regime up to the 
tricritical temperature. This barrier may be related to 
the interface tension for spin -0 droplets/strips, which we 
showed vanishes as one approaches the tricritical point 








from the first-order regime. Moreover, further numeri¬ 
cal simulations performed at the tricritical temperature 
Tt = 1.4182, proposed by Deserno 0 , indicated that 
this original estimate is rather accurate, verifying at the 
same time the expected Ising tricritical exponent values 
of V = 0.5 and ajv = \ from infinite-volume extrapola¬ 
tions of our effective exponents. 

A further asset of the proposed numerical and scal¬ 
ing schemes is that it opens a new window for revisit¬ 
ing the effect of disorder in first-order phase transitions 
in both two and three dimensions, where a unified ap¬ 
proach to universality is still missing. For instance, al¬ 
though it is known that in two dimensions under the 
presence of bond disorder the ex-second-order regime of 
the Blume-Capel model falls into the universality class of 
the corresponding random Ising model along the lines of 
the strong universality hypothesis [l^, the same is not 
true for the ex-first-order regime. Interestingly enough, 
for the ex-first-order regime different results have been 
obtained for different lattice geometries [3, [l^. The 
situation in three-dimensions in even more ambiguous, 
where one has to be also careful with respect to the dif¬ 
fused amount of disorder in the system in order to secure 
the switching to a continuous transition (d^ . |43| . A re¬ 
cent study of the random version of the three-dimensional 
Blume-Capel model suggested a possible new universality 
class at the ex-first-order regime, different to that at the 
ex-second-order regime fl9j |. an interesting finding if one 
considers that the two transitions are between the same 
ferromagnetic and paramagnetic phases. Yet, the authors 


of Ref. M clearly underlined the need for a more sophis¬ 
ticated approach (in both numerical and scaling terms) in 
order to tackle efficiently the low-temperature disorder- 
induced continuous transition regime of the model. 

To conclude, using as a platform the Blume-Capel 
model that shows the unique feature of having contin¬ 
uous and first-order transition lines in its phase dia¬ 
gram, we believe that the practise followed in the present 
manuscript applied over a wide range of disorder-strength 
values and temperatures will provide a better under¬ 
standing of the effect of disorder in spin systems. Us¬ 
ing the parallelized version of the multicanonical method 
and crossing the phase boundary along the crystal-field 
axis we expect to be able to study systematically the 
universality class and scaling corrections at the disorder- 
induced second-order phase transition of the Blume- 
Capel model, the shift behavior of the tricritical point 
as a function of the disorder strength, and other relevant 
open questions. Research in this direction is currently 
under way. 
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